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Coupled attitude and trajectory optimization for a launcher 

tilting maneuver 

Jiamin Zhu* Emmanuel Trelat^ Max Cerfl 


Abstract 

We study the minimum time control problem of the launchers. The optimal trajectories of 
the problem may contain singular arcs, and thus chattering arcs. The motion of the launcher is 
described by its attitude kinematics and dynamics and also by its trajectory dynamics. Based 
on the nature of the system, we implement an efficient indirect numerical method, combined 
with numerical continuation, to compute numerically the optimal solutions of the problem. 
Numerical experiments show the efficiency and the robustness of the proposed method. 


Keywords: Coupled attitude orbit problem; optimal control; Pontryagin maximum principle; 
shooting method; continuation; chattering arcs. 
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1 Introduction 

The usual way to approach the control of the guidance of satellites consists of considering 
separately the attitude and the trajectory dynamics. However, for the launchers like Ariane or 
Pegasus, the trajectories are controlled by their attitude angles, and it is desirable to determine the 
optimal control subject to the coupled dynamical system. Thus, we consider the time minimum 
control of the attitude reorientation of a rocket taking in consideration both the attitude movement 
and the trajectory dynamics. We denote this problem as (MTCP). Our objective is to design an 
efficient numerical strategy to solve the problem (MTCP). 

There are two main types of numerical approaches to solve an optimal control problem. On 
the one part, the direct methods (see, e.g., 0 ) consist of discretizing the state and the control and 
thus of reducing the problem to a nonlinear optimization problem (nonlinear programming) with 
constraints. On the other part, the indirect methods consist of numerically solving a boundary 
value problem obtained by applying Pontryagin maximum principle (PMP, see |10|1. by means of 
a shooting method. Both direct and indirect methods are not easy to initialize successfully, it is 
required to combine them with other theoretical or numerical approaches (see the survey m)- 
The numerical continuation is a powerful tool to be combined with the indirect shooting method 
based on the Pontryagin Maximum Principle. The idea of continuation is to solve a problem step 
by step from a simpler one by parameter deformation (see e.g. [1]). For example, in [H m H] , the 
continuation method is used to solve difficult orbit transfer problems. Here, we will use numerical 
continuation to numerically solve the problem. 

Indeed, the coupling of the attitude movement (fast) with the trajectory dynamics (slow) in 
the problem (MTCP) generates difficulties in numerical approaches, especially in the indirect 
methods, where the Newton-like methods are used to solve the boundary value problem. More¬ 
over, this property makes the numerical approaches extremely difficult to be initialized. Another 
difficulty in numerically solving the problem (MTCP) is the chattering phenomenon. Chattering 
means that the control switches an infinite number of times over a compact time interval. Such 
a phenomenon typically occurs when trying to connect bang arcs with a higher-order singular arc 
(see, e.g., [3 [Tills]). We proved in [Ml HI] that there exists such a chattering phenomenon in 
the problem (MTCP), and the chattering phenomenon is a bad news when applying numerical 
approaches. 

In view of these difficulties, we propose an efficient numerical continuation procedure to com¬ 
pute optimal and sub-optimal trajectories for the problem (MTCP). Due to the chattering phe¬ 
nomenon, numerical continuation combined with shooting cannot give an optimal solution to the 
problem for certain terminal conditions for which the optimal trajectory contains a singular arc of 
higher-order. In that case, our indirect approach generates sub-optimal solutions, by stopping the 
continuation procedure before its failure due to chattering. 

Our objective is to design and implement a numerical method having the following qualities: 
general (vehicle features, terminal conditions), robust, fast and automatic, so that it could be used 
efficiently by an inexperienced user on a large scope of applications of the real world. Numerical 
experiments indicates that this approach happens to meet to all the expected features. 

The paper is organized as follows. In Section [3 we establish the model of the rocket move¬ 
ment and we formulate the problem (MTCP). In Section [3 the PMP is applied to the problem 
(MTCP) and some theoretical results on the extremals are introduced. In Section[4l two auxiliary 
problems are introduced and the continuation procedure is stated. Some numerical tips improving 
the proposed numerical strategy are also presented. Finally in Section [3 numerical results are 
given. 


2 


2 Problem Statement 


We will consider two types of launchers: one is a classical launcher like Ariane rocket, the other 
one is an airborne launcher like Pegasus rocket (being launched horizontally from an airplane). 
In the case of an Ariane-type launcher, the aerodynamical forces are very small compared with 
its thrust and the gravity, and so we ignore them. However, in case of a Pegasus-type launcher, 
we will need to add the aerodynamical forces to the model. Throughout the paper, we make the 
following assumptions: 

• The Earth is a sphere and is fixed in the inertial space. 

• The position and the mass of the rocket are constant during the maneuver. 

• The rocket is an axial symmetric cylinder. 

• The rocket engine cannot be shut off during the maneuver and the module of the thrust force 
is constant, taking its maximum value, i.e., T = T^ax- 

2.1 Coordinates and Model 

All the coordinate systems introduced here are Cartesian coordinate systems. 

The launch frame (reference frame) Sr = {xR,yR, zr) is fixed around the launch point Or 
(where the rocket is launched). The axis xr is pointing radially outwardly (normal to the local 
tangent plane), and the axis zr points to the North. 

The body frame Sb = {xb, yb, Zb) is defined as follows. The origin of the frame Ob is fixed 
around the mass center of the rocket, the axis Zb is along the axis-symmetric axis of the rocket, and 
the axis Xb is in the cross-section. The body frame can be derived by three ordered unit single-axis 
rotations from the launch frame, 
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where 9 is the pitch angle, ip is the yaw angle, (p is the roll angle and Ra{b) means to rotate the 
frame around the axis a G {x, y, z} by an angle 6 G M. Therefore, the transfer matrix from Sr to 
Sb is LbR = R4(P)Rxiip)Ry{e). 

The velocity frame Sy = {xv,yv, Zy) is fixed around the mass center of rocket. The axis Xy 
is parallel to the velocity of the rocket v, and the axis Zy is normal to the velocity, pointing to 
the direction of the lift force L of the rocket. This frame can be derived by two unit single-axis 
rotations from the launch frame, 


Sr 


Rx{^) Cl 

-O -S’ b, 


where ^ is the flight path angle and k is the bank angle of the flight. Denote by v the module of 
the velocity vector, i.e., v = ■sjvP;. + Vy + v^. since {v)r = (vcos^, usin^sin/t, —vsin^cos k)^, we 
have 

COS^ = Vx/v, tSUlK = —Vy/Vz- 

This frame will be used to introduce aerodynamical forces (lift L and drag D) into the model. 


Attitude motion The dynamical and kinematic equations of the attitude are 

9 = {ujx sin (p + ujy cos </>) / cosip, ip = uix cos (p — ujy sin (p, p = (ujx sin p + ujy cos p) tan ip, 

UJx = —bU 2 , dly = bui- 


( 1 ) 
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where u = {ui,U 2 ) G is the control input expressed in frame Sb, (cJ)b = (wxjUJy) is the angular 
velocity of the launcher cJ expressed in frame Sb (assuming that there is no rotation along the 
asymptotical axis of the launcher), and b is the angular acceleration constant depending on fimax 
(the maximum angle between the thrust T and the axis Zb), Ic (the inertia around the axis Xb), I 
(the length of the rocket). Note that ||it|| = ^ 1. See m for details of this model. 


Trajectory dynamics The drag D and lift L are calculated by = —^pv^SCxXy, L = 
^pv'^SCxZy, where p is the air density, S is the reference surface, Cx and are constant co¬ 
efficients of the drag and the lift. In the velocity frame Sy, the components of the vectors D and 
L are expressed as {D)v = {—Cx'mv'^,0,0)^ , {L)y = (0,0, where Cx = ^pSCx/m and 

Cz = \pSCx/m. By using the transfer matrix from Sy to Sr, we get 

{D)r = Cxmv^{— cos^, — sin^ sin k, sin^ cos k)^, {L)r = Czmv^{— sin^, — cos^ sin k, cos^ cos k)"'', 
and so we get the trajectory dynamical equations in frame Sr as 
Vx = a sin 6 cos ij} + px — CxV^ cos ^ sin 

Vy = —asintf} + Py — CxV^ sin^ sin«; — CzV^ cos^ sinK, (2) 

Vz = a cos 0 cos ijj + Pz + CxV^ sin^ cosk -|- CzV^ cos^ cosk, 

) is the velocity of the launcher v expressed in Sr, {p)r = {px,Py,Pz) is 


'Xl ^y, 


where (v)r = {t 

the gravity vector expressed in Sr, a = Tmax/m is the thrust acceleration. For the Ariane rocket, 
we do not consider the aerodynamical forces, i.e., Cx = Cz = 0. 


2.2 Minimum Time Control Problem (MTCP) 

Defining the state variable x and the control variable u as 

X = {Vx,Vy,Vz, 0 ,'llj,(l),UJx,UJy), U = {ui,U 2 ). 

The initial conditions are defined by Vxq , Vy^, Vz^, 0o, V’o: </'o; and cuy^, 


= ( 0 ) = 


/(0)=' 


'yo^ 




^(0)=bo, V'(0)=V'o, (/>(0) = (/> 0 , ^^x( 0 )=UJxo, UJy{ 0 )=UJy 


( 3 ) 


The desired final velocity is required to be parallel to the axis Zb, according to (F(t/))_R A 
{zb{tf))R — 0. The constraints on the final conditions are defined by Of, ipf, (pf, uixj and uiyj. 


Vzf sin Ipf + Vyj cos Of cos'0/ = 0, Vzf sin Of — Vxf cos Of = 0, 

e{tf)=0f, P){tf)='iPf, (P{tf)=(Pf, ujx{tf) =u;xf, ujy{tf)=i 


( 4 ) 


^Vf 


Note that the parallel condition on the final velocity is due to the fact that most rockets are planned 
to maintain a zero angle of attack along the flight. The angle of attack, when the wind is set to 
zero, is defined as the angle between the velocity and the rocket body axis. 

We set Xq = {vxo,Vy^,VzQ,Oo,ipo,(po,u!xo,ujyg), and we define the target set 

Ml ={{vx,Vy,Vz,0,'ip,(p,uJx,u}y) I VzSinipf + Vy cosOf costpf = 0, 

VzSinipf + VyCOsOf costpf = 0, 0 = 0f, 'tp = ipf, cp = (pf, ojx = uJxf, ujy=ujy^}. 

The minimum time control problem (MTCP) consists of steering the system (II|)-(l2]) from a;(0) = xq 
to the final target Mi in minimum time tf, with controls satisfying the constraint uf + U 2 ^ 1. 

For convenience, the velocity will be defined by polar coordinates, i.e., Vx = v sinOy cosipy, 
Vy = —vsinipy, Vz = V cosOy cosipy, where Oy and ipy are “pitch” and “yaw” angles of the velocity 
vector. 
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3 Theoretical Results 

According to the Pontryagin Maximum Principle (PMP), the Hamiltonian of the optimal control 
problem (MTCP) is 


H{x{t),p{t),p°,u{t)) = {p{t),f{x{t),u{t))) +p°, 

where p = {pvxTPvy,Pvz,P 0 ,Piji,P<t,,Puix,Puy) is the adjoint variable, f{x,u) is the derivative of the 
state variable x, i.e., x = f{x,u) and ^ 0 is a real number. Here we assume = —1 (see 
m)- The differential equation of the adjoint variable p is given by the partial derivative of the 
Hamiltonian as 

p{t) =- — {x{t),p{t),p°,u{t)). (5) 

The maximization condition of the PMP yields, almost everywhere on [0,t/], 

^h,{ty + h2{ty ll$(0r 

whenever = {hi{t),h 2 {t)) ^ (0,0). Here we have hi{t) = bpcuyit) and h 2 {t) == —bpujx{t)- We 
call $ (as well as its components) the switching function. 

Moreover, we have the transversality condition p(t^) T where is the tangent 

space to Ml at the point x{tf), i.e., 

Py^ sini/jf = Pv^ sinOfCOSipf + Pv, cos 9 f cos tpf, ( 6 ) 

and, the final time tf being free and the system being autonomous, we have also ho{x{t),p{t)) + 

||$(t)|| +p° = 0 , v< e [o,t/]. 

The quadruple {x{-),p{-),p^,u{-)) is called an extremal lift of a;(-). An extremal is said to be 
normal (resp., abnormal) if < 0 (resp., = 0). We say that an arc (restriction of an extremal 
to a subinterval I) is regular if ||$(t)|| ^ 0 along I. Otherwise, the arc is said to be singular. 

A switching time is a time t at which $(t) = (0,0), that is, both hi and h 2 vanish at time t. 
An arc that is a concatenation of an infinite number of regular arcs over a compact time interval 
is said to be chattering. The chattering arc is associated with a chattering control that switches 
an infinite number of times, over a compact time interval. A junction between a regular arc and a 
singular arc is said to be a singular junction. 

When Ca; = 0 and Cz = 0, the explicit form of equation ([S]) is given by 


Pv^ = 0, Pvy = 0, Pv^ = 0, 

P 9 = —acosijj{py^ cos9 — py^ sin6*), 

PjP = a sin Ip sin 9 py^ + acosippy^ + acos9sintppy^ — sintp{u}x sini/) + uiy cos(p)/ cos^ ippg 

- {ojxfii'o-cp + ujycoscp)/cos^ ipp^, 

= — {ujx cos (p — LOy sin (p) / cos ippe + {ujx sin (p + ujy cos (p) 

— tan'ip(ujx cos (p — ujy sin (p) p^, 

Puj^ = —sin(p/ cosippe — coscpp^jj — sin'i/’sin^/cos'i/’p^, 

Pujy = — coscp/ cosippe + sin(pp.li, — sincos <p/ cosipp^j,. 

We have shown in m that the singular extremals of the problem (MTCP) are normal ones, i.e., 
p° ^ 0, and that they are of intrinsic order two. In that case the singular junction can only be 
realized by chattering, i.e., if the singular junction happens at time r and the control is regular 
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on (ti,T) and is singular on (r,^ 2 ), then the control has to switch an infinite number of times on 
(r — e,r), Ve G K such that 0 < e < — r. 

The concatenation of regular arcs are classified in m by their contact with the switching surface 
(filled by switching points), and we know that if the switching points of the problem (MTCP) are 
of order one, then the control will turn an angle tt when passing the switching surface. We will see 
this phenomenon in the numerical results. 

Actually, when ^ 0 and ^ 0, the derivative of = [pvx^Pvy^Pvz) is no longer zero. 
However, the terms induced by the air draft and lift do not change the Lie bracket conhguration 
of the system, and thus the results obtained in [15] are still valid. 

Note again that the chattering causes hard problems when applying numerical methods. In that 
case, we can no longer obtain an optimal solution numerically. However, a sub-optimal solution 
can be derived by stopping the continuation procedure proposed in Section 14.21 before chattering 
occurs. 


4 Auxiliary Problems and Continuation Procedure 

To construct a suitable continuation procedure, we introduce two auxiliary problems. The 
first one is the problem of order zero, denoted by (OCPO). The solution of this problem can 
be easily computed. Then we just need to plug this simple, low-dimensional solution in higher 
dimension, in order to initialize an indirect method for the more complicated problem (MTCP). 
The second one is the regularized problem, which is used to overcome the chattering issues. Since the 
solution of problem (OCPO) is contained in the singular surface (filled by the singular solutions) 
for the problem (MTCP), passing directly from the solution of problem (OCPO) to the problem 
(MTCP) makes the optimal extremals to contain a singular arc (and thus chattering arcs), and the 
shooting method will certainly fail due to the numerical integration of discontinuous Hamiltonian 
system. 


4.1 Two auxiliary problems 

Problem of order zero. We define the problem of order zero, as a “subproblem” of the complete 
problem (MTCP), in the sense that we consider only the trajectory dynamics (without aerody¬ 
namical forces) and that we assume a perfect control meaning that the attitude angles (Euler 
angles) take the desired values instantaneously. Thus, the attitude angles are considered as control 
inputs in that simpler problem. Denoting the rocket axial symmetric axis as e and considering it 
as the control vector (which is consistent with the attitude angles 9, ■0), we formulate the problem 
as follows: 

v = ae + g, v{0) = vq, v{tf)//w, ||u;|| = 1 , mint/, 

where w is a given vector that refers to the desired target velocity direction. This problem is easy 
to solve, and the solution derived by applying the PMP is 


1 / kw — Vo 

a V tf 


-02 -I- - 4ai03 

=- 2 ^^-’ = 




with k = {vo,w) + {g,w)tf, ai = of - \\{g,w)w - gW^, 02 = 2{{vo,w){g,w) - (uo,g)), and 03 = 
-||{uo,u;)u;-uoll^. 

Since the vector e is expressed in the launch frame Sr as {e)R = (sin 9 cos ip, — sin 0, cos 9 cos 0)^, 
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the Euler angles 9* G (—7r,7r) can be calculated by assuming tp* G (—7r/2,7r/2) as 


( arctan(|e*/e*|), e* > 0, 

6* = atan2(e;, e*) = sign(e;) 7r/2, e* = 0, (8) 

[(tt - arctan(|e;/ej|)), e* < 0, 


where sign(-) is the sign function, and 


tp* = arcsin(—e^. (9) 

The assumption of tp* G (— 7 r/ 2 , 7 r/ 2 ) generally works for the Ariane-type rockets since the ma¬ 
neuvers are generally quasi planar and the change in the yaw angle is small. However, for the 
Pegasus-type rockets, which could make large yaw angle maneuvers, we have to consider also the 
case cos tp* < 0 leading to 

0 * = atan 2 (—e*, —e*), i/;* = —sign(ep( 7 r — arcsin(| — e*|)). ( 10 ) 

Both [TUI and ([5]) and ([HI) are possible solutions. We will choose the pair {0* ,tp*) minimizing the 
value |' 0 o — ' 0*1 + l' 0 / ~ ' 0 * 1 - 

Indeed, when tp = ± 11/2 + kir, fc G Z, the Euler angles are not well defined, and if tpo and tpf 
are given near ± 7 r/ 2 , we need to “change” them to a be near ± 7 r/ 2 . This is done in section [4?3l 
below by transforming of the reference frame Sr to a new reference frame The singularities 
of Euler angles are also treated in section 14.41 by calculating the limit values of the vector field at 
these singular points. 

Given a real number 0*, the optimal solution of the problem (OCPO) actually corresponds to 
a singular solution of (MTCP) with the terminal conditions given by 

Ux(0) = Vy{0) = Vyg, v^{0) = 

0 ( 0 ) = r, 0 ( 0 ) =00 , 0 ( 0 ) =00 a;,( 0 ) = 0 , a;y( 0 ) = 0 , ^ ^ 

Uz(t/) sin0/ -I- Uy(t/) COS0/COS0/ = 0, Vz{tf)siTt9f—Vx{tf)cos6f = Q, (12) 

0{tf) = 0*, tp{tf)=tp*, ,cP{tf)=(P*, ujx{tf) = 0, ujy{tf)=0. (13) 

It is worth nothing that this solution is lying in the singular surface of the problem (MTCP) 
meaning that it is on the “highway” between two fixed points in the state space. On this “highway”, 
the system goes the most rapidly towards to aimed final point or manifold. Indeed, we observe 
in the numerical simulations that the singular extremals of the problem (MTCP) acts a role 
similar to that of the stable points in the “turnpike” phenomenon as described in m- the optimal 
trajectories first tend to reach the singular surface (to have a greater speed in transferring its state), 
then stay in the singular surface for a while until it is sufficiently close to the target submanifold 
Ml, and finally get off the singular surface to join the target submanifold. Note that the singular 
arc is not necessary if the state is sufficiently close to the target. 

In view of (HUl-IIll), a natural continuation strategy is to simply vary the terminal condi¬ 
tions step by step until they correspond to Mi. However, as we have mentioned, the chattering 
phenomenon causes the failure of this simple strategy, and thus we introduce another auxiliary 
problem, the regularized problem, in which the singular arcs of the problem (MTCP) become 
regular arcs, and thus the difficulty caused by chattering is bypassed. 
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Regularized problem. Let 7 > 0 be arbitrary. The regularized problem (OCPR)..^ consists of 
minimizing the cost functional 

ftf 

Cj=tf+j {ul + U2)dt, (14) 

Jo 

for the system under the control constraints —with terminal conditions 

©- 0 . Note that, here, we replace the constraint uf + uf ^ 1 with the constraints — 1 ^ Ui ^ 1. 
The advantage, for this intermediate optimal control problem with the cost (HI, is that the 
extremal controls are continuous. 

The Hamiltonian is 

H-y = (p,fix,u))+p°{l+-ful+jul), (15) 

and according to the PMP, the optimal controls are 

ui{t) =sat(-l,- 6 p^^(t)/( 27 p°),l), U 2 it) = sat{-l,bp^^{t)/{2jp°),l), (16) 

where the saturation operator sat is defined by sat(— 1 ,/(t), 1 ) = —1 if f{t) ^ — 1 ; 1 if f{t) ^ 1 ; 
and fit) if -1 ^ f{t) < 1 . 

Indeed, an extremal of (OCPO) can be embedded into the problem (OCPR).y by setting 

M(t) = (0,0), e{t)=9*, = (j)it)=4>*, ujj:it)=0, ujy{t)=0, 

peit) = 0, p^it) = 0, p^(t) = 0, Pu:xit) = 0, Pa;yit) = 0, 

where 6* and tf* are given by dl]) and o, with terminal conditions given by (HB-dll. Moreover, 
it is not a singular extremal for the problem (OCPR).y . The extremal equations for (OCPR).y 
are the same than for (MTCP), as well as the transversality conditions. 

Note that the difficulty caused by the chattering still exists when trying to go from the regular¬ 
ized problem back to the problem (MTCP) if there exists a singular arc in the optimal trajectory. 
However, in that case, by stopping at some point during the last continuation step (see below), an 
acceptable sub-optimal trajectory will be found for the problem (MTCP). 

4.2 Continuation procedure 

The ultimate objective is to compute the optimal solution of the problem (MTCP), starting 
from the explicit solution of (OCPO). 

Numerical strategy We proceed as follows: 

• First, we embed the solution of (OCPO) into (OCPR).y . For convenience, we still denote 
by (OCPO) the problem (OCPO) seen in high dimension. 

• Then, we pass from (OCPO) to (MTCP) by means of a numerical continuation procedure, 
involving four continuation parameters: the first three parameters Ai, A 2 and/or A 3 are used 
to pass continuously from the optimal solution of (OCPO) to the optimal solution of the 
regularized problem (OCPR).y , for some fixed 7 > 0, and the last parameter A 4 is then 
used to pass to the optimal solution of (MTCP) (see Figure [T]). 

The unknowns of the shooting problem are Pv^iO), PvyiO), Pv^iO), Pe(0), p^iO), ^^(0), Pujy:i0), 
Pujy (0) and tf . When D = 0 and L = 0, we have that py ^, py^ and py^ are constants, hence by using 
the transversality condition ([ 6 ]) , we can easily reduce the number of the unknowns by one. In view 
of more general cases, we will not do this reduction and we will use the transversality condition as 
a part of the shooting function. In the following, we denote the continuation step corresponding 
to parameter A^ by Ai-continuation. 



Figure 1: Continuation procedure. 


Ai-continuation and A 2 -continuation The parameter Ai is used to act, by continuation, on 
the initial conditions, according to 

0(0) = r(l-Ai) + 0oAi, V^(O) =?/'*(l-Ai) + V'oAi, 0(0) = <^*(l-Ai) + .^oAi, 

<j-’a;(0) = a;*(l — Ai) + Waj^Ai, Wy(0) = a;*(l — Ai) + uiy^Xi, 

where w* = w* = 0, 0* = 0, and 0*, 0* are calculated through equation ([S])-®. 

The shooting function S\^ for the Ai-continuation is defined by 

5 'ai = (pujAtf), Pe{tf), P^(.tf), p^{tf), 

u^(t/) sin0/ + Vy{tf) COS0/ cos0/, Vz{tf) sin0/ — Vx{tf) cos0/ 
sintpf — {py^ sm9f cos0/ + Pv^ cos0/ cos0/)), 

where with = — 1 is calculated from (ITSl) and Ui and are given by (|TO1) . 

Note that we can use S\^ as shooting function owing to (OCPR)..y . For the problem (MTCP), 
if S'ai = 0 , then together with ojxitf) = 0 and ojy{tf) = 0 , the final point {x(tf),p{tf)) of the 
extremal is then lying on the singular surface and this will cause the failure of the shooting. 
However, for problem (OCPR).^ , even when x{tf) belongs to the singular surface, the shooting 
problem can still be solved. 

Initializing with the solution of (OCPO), we can solve this shooting problem with Ai = 0, and 
we get a solution of (OCPR).^ with the terminal conditions (ITT]) - (fT^ (the other states at tf being 
free). Then, by continuation, we make Ai vary from 0 to 1, and in this way we get the solution 
of (OCPR).y for Ai = 1. With this solution, we can integrate extremal equations o, © and 0 
to get the values of the state variable at tf. We denote 0e := d{tf), 0e := fpitf), 0e := 0(^/), 
ujxe '■= <-^x{tf) and lUye '■= (^y{tf) the “natural” conditions obtained at the final time. 

In a second step, we use the continuation parameter A 2 to act on the final conditions, in order 
to make them pass from the “natural” values 0e, tpe, 0e, tOxe and coye, to the desired target values 
0/, tpf, (j)f, LOxf and uiyf. The shooting function is 

S\2 — (uJx{tf^ (1 X2')UJxe X2UJxfj (1 X2^iOy^ A2CUy^, 

9{tf) - (1 - \2)9e - A 26 */, lj}(tf) - (1 - A2)0e “ A20/, 0(t/) “ (1 “ A2)0e - A20/, 
Vz{tf) sin0/ 'Vy{tf) cos 0/ cos 0/, Vz{tf) sin0/ — Vx{tf) cos0/, 

Py^ sin0/ - {py^ sin0/ COS0/ + Py^ cos0/ cos0/, H^{tf)). 

Solving this problem by making vary A 2 from 0 to 1, we obtain the solution of (OCPR)-j, with the 
terminal conditions ®- 0 - 

As-continuation The parameter A 3 is added to the terms induced by the aerodynamic forces 
in system 0 - 0 , i.e., replace Cx and Cz by X^CxD and AsC^. It is the only parameter that acts 
on the differential system, and it is used to differ the Ariane-type and Pegasus-type rockets. We 
let A 3 = 0 during the Ai-continuation and the A 2 -continuation, which means during these first 
two steps of the continuation procedure, we do not consider the influence of the drag and the lift. 
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Then, if we are dealing with a Pegasus-type rocket, we vary A3 from 0 to 1 to derive a solution of 
problem (OCPR)-,, with drag and lift, and the shooting function remains the same as for the A2- 
continuation. Otherwise, for an Ariane-type rocket, we skip Aa-continuation, i.e., we keep A3 = 0 
and go ahead to the A4-continuation. By doing this, the program is suitable for both types of 
rockets. 


A4-continuation Finally, in order to compute the solution of (MTCP), we use the continuation 
parameter A4 to pass from (OCPR).^ to (MTCP). We add the parameter A4 to the Hamiltonian 
and to the cost functional (iMll as follows: 

j-tf 

Cj=tf+^ {ul + ul){l - X 4 ) dt, 

Jo 

H{tf,X 4 ) = {pj{x,u)) +p° +p°-f{ul+ul){l - X 4 ). 

Then, according to the PMP, the extremal controls are given by Ui = sat(— 1 ,Uie, 1 ), i = 1 , 2 , 
where 

bPcjy -bpuj^ 

-2^07(1 - A4) -f hX 4 4 jpl^ +pI^ -2p0^(l - A4) -h hX 4 4 jpl^ +pI^ 

The shooting function Sx^ is the same as replacing Hj{tf) with H^{tf,X 4 ). The solution of 
(MTCP) is then obtained by making vary A4 continuously from 0 to 1 . 

Note again that the above continuation procedure fails in case of chattering, and thus cannot 
be successful for any possible choice of terminal conditions. In particular, if chattering occurs then 
the A4-continuation is expected to fail for some value A4 = A| < 1 . But in that case, with this 
value A4, we have generated a sub-optimal solution of the problem (MTCP), which appears to 
be acceptable and very interesting in practice. Moreover, the overall procedure is very fast and 
accurate. Note that the resulting sub-optimal control is continuous. 

As we have mentioned, it is possible that an optimal trajectory meets the singularities of the 
Euler angles, i.e., cosijj{t) = 0 , for some t S [ 0 ,t/]. A coordinate system transformation can be 
made to change the set terminal conditions and avoid meeting the singularity. We will use two 
numerical tricks as follows. The first one is to use a new reference frame S'^, instead of Sr, such 
that the terminal conditions in the new reference frame are easier to solve. The new frame is 
set by rotating the initial frame Sr. This amounts to a nonlinear state transformation, and leads 
to a preconditioner that makes the proposed continuation procedure more robust. The second trick 
is to overcome the singularities of the Euler angles by redefining the vector fields at the singular 
points. 


4.3 Change of frame 


We define a new coordinate Sr' which is derived by two single-axis rotations from the frame 
Sr, i.e., 

Ry{a) iia,(/3) 

Sr -S- o -S- Sri, 


and so the transfer matrix from Sr to S'^ is Lrir = Rx{f3)Rz{a). Denoting the Euler angles of Sb 
with respect to the new frame Sri as 9', ip' and (j)', we get the transfer matrix from Sf^ to Sb as 
LbRi = RziP')RA^’)Ryi9’). 

Therefore, it follows from LbRLRRi — Rz{4>)Rx{'tp)Ry{9)L^, = LbRi that the angles 9', ip' and 
4>' are functions of the angles 9, ip, (p, a, and (3. Moreover, the velocity v in the S'j^ can also be 
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obtained by (■y)^ = Lii{v)n. The angular velocity uj = (uJx,ojy) is expressed in the body frame 
Sb and it is therefore not altered by this coordinate change. 

Given fixed a and (3, the change of frame corresponds to a nonlinear invertible change on the 
state variable as follows 

x' = 

where ifatti’) is the mapping from the Euler angles in Sr to the Euler angles in and Id is 
identity matrix of order 3. Note that this transformation is invertible. For every given value of a, 
the value of /3 can be chosen such that '0^(t/) = ij^'f = ~V'o = By doing this, the terminal 

values on ip is “closer” to the origin and hence “farther” from the singularities of the Euler angles. 
Indeed, this only reduces the risk of meeting with the singularities of the Euler angles, so we will 
as well use this numerical trick in the next section to improve the robustness of the numerical 
procedure. 

Further, we explain why this state transformation can be seen as a preconditioner to our 
numerical method. In view of the shooting method, we use a Newton-like method to solve the 
boundary value problem. For simplicity of explanation, we use the simplest Newton method and 
denote the variables of the shooting method as z. Given an initial guess oi z = Zq, the equation 
S'(z) = 0 is solved iteratively, i.e., 

Jizk)zk +1 = J{zk)zk - S'(zfc). 

where J = dS/dz. Define a diffeomorphism ip{-) such that z = Then the original problem 

becomes solving S{y) =0 and the Newton method iterative steps become 

J{yk)yk+i = J{yk)yk - S{yk), 


where 

Qz 

J{yk) = J{zk)^{yk)- 

It is easy to see that the matrix M = bere actually acts as a preconditioner in the shooting 

method and it can be used to get a Jacobian matrix J with a smaller condition number. In the 
numerical experiments, we use a Fortran subroutine hybrd.f (see [5]) which uses a modification of 
the Powell hybrid method: the choice of the correction is a convex combination of the Newton and 
scaled gradient directions, and the updating of the Jacobian by the rank-I method of Broyden. 

Note in addition that the differential equations for the new variable x' keep the same form as 
the old variable a;, and by using the PMP, the adjoint vector p' to the new state x' can also be 
derived from {x,p), i.e.. 


/ T t ' I ' \X /dpatt(x)— 1 

Pv = Lr'RPv, [P9,P,p,P^) = ( -- ) 


(Pe.Pij.P^V, (pL.pLJ = {Pu::,,Pi^y), 


and this is also invertible. 

To sum up, the reference frame can be chosen such that the problem (MTCP) can be easier 
solved numerically. However, a priori, we do not know which pair of (a, /3) is the most suitable 
to the problem. If we take the pair that makes the condition number of the Jacobian matrix 
J at the initial time smallest, the value of |0g — 9'j\ may become too large (huge pitch angle 
maneuver) leading to the fail of the numerical strategy and it does not ensure that the condition 
number remains small in all the four continuation steps. If we choose the pair that makes the 
value \9q — the smallest, it is possible that the Jacobian matrix is ill-conditioned. Therefore, 
we propose to do the following: 

Frame change heuristics 
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• stepl: fix arbitrary a, calculate /3 so that ipf = ~' 4 ’o 3 'iid assess the new terminal conditions; 

• step2: solve the problem (MTCP) using the continuation procedure proposed; 

• step3: if the step2 succeeds, stop. If the step2 fails, repeat stepl-step2 with a new a = a + 5a 
{Sa constant). 

We will see in the numerical results that this enhances the robustness of our numerical methods 
though it consumes more time. 

4.4 Singularities of the Euler angles 

As we can see from o, when Ip = TV 12 + kn, fc S Z, the system encounters singularities due 
to the fact that the Euler angle (p is not well defined at these points. Here we call these points 
singularities of the Euler angles. One way to overcome these singularities is to calculate the limit 
values of the differential of x and p at these points. Assume that 6 is bounded, then we have 
{lOx sin (p + Uy cos (/))—>■ 0 when i/) —^ 7 r /2 + kn. Thus we get from 

lim 9p = lim {ujx sin cp + ujy cos sin ijj = 0, lim OI(f>=l, 

'il)—¥7Tj2-\-kn i/j—)-7r/2+fc7r j 2 -\-k'K 

that 9 = <j> = 0 when ip —^ 7 r /2 + kn. Assume that lim,^^^/ 2 +fe 7 r = A < oo, then it 

follows from 

A — lim Pe + sin ip _ pg + sin ip + cos iptp _ 

/2+kTT COS'tp i/j—>7r/2+fc7r sinpj'tp 

that A = 0. Hence, we get 

Pe = 0 , P 0 = 0 , p^ = a sin 9pyx + a cos 9pyz , Pujx = -Pi/- cos (p, Pu,y = p^ sin p. 
Summing up, at points ip —)• tt /2 + kn, equation ([T]) become 

9 = 0, Ip = LOx cos p — LUy sin p, p = 0, Cox =—hu 2 , ujy = bui, (17) 

and 0 become 

P 8 = 0, p^ = asin9pyx + acos0p„z, p^ = 0, Puix = -Pip cosp, p^y = p^ sinc^i, (18) 

When we are close to such a singularity, we use (HU and (UHl) instead of the original system 
equations O- 


5 Numerical results 

5.1 Statistical tests 

We test the proposed numerical method in three typical cases of launchers in order to check its 
robustness with respect to the terminal conditions. The first case is the atmospheric ascent phase 
of the Ariane-type rocket, the second one is the flight phase of the Ariane-type rocket (after the 
first stage separation), and the third one is the launch phase of the Pegasus-type airborne rocket. 
The launcher data (see equations 0-@) in these three cases are given in the Table [H We also 
indicate the computer features on which the tests were run. 
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Ariane launch 

Ariane flight 

Pegasus 

System parameters 

a = 18, 6 = 0.0138 

a = 25, 6 = 0.0158 

a = 24.873, 6 = 0.0607 

Aerodynamical parameters 

Cx — C 2 : = 0 

Cx — Cz — 0 

Cx = Cz = 5 X 10“'^ 

Machine info 

CPU: Intel Core i7-4790S CPU 3.2GHz Memory: 7.6 Gio 
Compiler: gcc version 4.8.4 (Ubuntu 14.04 LTS) 


Table 1: Launcher data. 

We solve the problem (MTCP) with different terminal conditions. The initial and final condi¬ 
tions are swept in the range given in Tabled The two last lines of the table define the restrictions 
applied to the terminal conditions in order to exclude unrealistic cases. 



Ariane launch 

Ariane flight 

Pegasus 

^^0 

[50, 500] m/s 

[2000,6000] m/s 

[300,300] m/s 

dvO 

fixed 90° 

[0,60]° 

[- 10 , 0 ]° 

i’vO 

fixed 0 ° 

fixed 0 ° 

fixed 0 ° 

00 

fixed 90° 

[0,60]° 

[0,50]° 

■00 

fixed 0 ° 

fixed 0 ° 

fixed 0 ° 

0f 

[60,85]° 

[0,60]° 

[60,90]° 

ijjf 

fixed 0 ° 

fixed 0 ° 

[0,90]° 

^xO 

fixed 0 ° 

fixed 0 ° 

[- 10 , 10 ]° 

OJyO 

fixed 0 ° 

fixed 0 ° 

fixed 0 ° 

0 

1 

free 

[- 20 , 20 ]° 

free 

0vO — 00 

free 

fixed 0 ° 

[- 10 , 0 ]° 


Table 2: Parameter ranges. 

For each variable, we choose a discretization step and solve all the possible combinations of this 
discretization (factorial experiment). There will be respectively 108, 234, 480 cases for the three 
launcher configurations. 

Results without change of frame The statistical results without change of frame are given 
in Table 131 We see that the success rate for the Ariane launch case (85.4%) is much higher than 
the Ariane flight case (50 %) and the Pegasus case (24.6 %). 

The table presents the average time spent for each continuation stage (Ai to A 4 ) and the average 
number of simulations performed to achieve each stage. These average performances are assessed 
on the successful cases. 

The main reason for this phenomenon is the singular arc (which is considered as the “highway” 
in our problem), which causes chattering. This can be seen in the following aspects: 

• The first one is the ratio of parameters a and b. By comparing the case of Airane launch 
and the case of Pegasus, we find that when a /6 is large (Ariane launch), the continuation 
procedure works better. For a launcher, the parameter a represents the capacity to reorient its 
velocity, whereas the parameter 6 represents the capacity to reorient its attitude. Hence, when 
a /6 is small (Pegasus), it means that the attitude can be turned quickly to the “highway”, 
which is the singular surface, and that there will probably be a singular arc in the optimal 
trajectory. 

• The second one is the velocity. In case of Ariane flight, though its ratio a /6 is large (even a 
bit larger than Ariane launch case), the velocity is much bigger than in the Ariane launch 
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case. Reorienting the velocity by the same angle requires more time than for a small velocity 
(see M)- That is to say, in order to reorient the velocity vector to the right direction, in 
case of Ariane flight, it needs to stay more time in the “highway”. 

• The third one is the reorientation angle. In the case of Pegasus, though the velocity is small, 
the difference between the initial Euler angles and the final Euler angles are large compared 
with the two other cases. This means that the velocity also has to have a large change 
in direction, which leads again the optimal trajectory to have high possibility to contain a 
singular arc. 


Statistical results without change of frame 


Ariane launch 

Ariane flight 

Pegasus 

Number of cases 

108 

234 

480 

Rate of success (%) 

- Total 

82.4 

50.0 

24.6 

- Before A 4 -continuation 

82.4 

75.6 

40.2 

Average execution time (s) 

- Total 

3.44 

25.66 

26.25 

- In Ai-continuation 

0.16 

2.52 

0.60 

- In A 2 -continuation 

2.63 

6.30 

3.35 

- In As-continuation 

0 

0 

2.46 

- In A 4 -continuation 

0.62 

18.23 

14.54 

Average number of simulations 
- Total 

122 

195 

188 

- In Ai-continuation 

20 

37 

47 

- In A 2 -continuation 

83 

60 

37 

- In Aa-continuation 

0 

0 

30 

- In A 4 -continuation 

19 

98 

74 


Table 3: Statistical result without change of frame. 

In term of execution time, it appears that most time is consumed in the A 4 -continuation (see 
Ariane flight and Pegasus cases). So far, we say that the method is efficient for the Ariane launch 
case, but it is less efficient and robust for the two other cases. Now we run the same tests using 
the change of frame. 

Results with change of frame The statistical result with change of frame are presented in 
Table m The change of frame is defined by the angles (a,/3). The resolution is restarted with 
varying (a,/3) values until successful Isee 14.31 Frame change heuristics). The pair {a*, (3*) is the 
pair that makes the continuation procedure works. In Table |4l the total execution time includes 
the {a, j3) search, while the continuation times account only for the final resolution once (a*,/3*) 
have been found. 

It is obvious that the proposed continuation procedure is more robust when combined with 
change of frame. Especially in the Pegasus case, the total success rate is improved from 24.6% to 
65%. 

Again, the success rates of Ariane flight and of Pegasus are lower than those of the Ariane 
launch, and the A 4 -continuation consumes more time than the three other continuations. 

In fact, when the optimal trajectory does not contain chattering arcs, the total execution time 
is just about several seconds, however, when the chattering phenomenon tends to appear (by 


14 











varying A 4 from 0 to 1 ), the continuation procedure start to take smaller and smaller steps and 
the execution time becomes much longer. This is why A 4 -continuation takes more time. 

Therefore, if we stop the A 4 -continuation when chattering occurs, we can obtain a sub-optimal 
solution within a very short time. Moreover, the success rate for obtaining a sub-optimal so¬ 
lution (success rate without A 4 -continuation) is much higher (88.9% for Ariane launch, 88 % for 
Ariane flight and 83.5% for Pegasus). This indicates again that most failures happen in the A 4 - 
continuation. 

We observe that either the solution is obtained very quickly or it can take a long time. Therefore, 
for a statistical analyses, we define arbitrarily two group of cases, denoted “easy” and “difficult”. 
The easy cases are those solved in less than 50 s, whereas the difficult cases are the other ones. 

We see that once the right reference frame is chosen, the time for solving the problem is 
reasonable even in the difficult cases. Much time is consumed by trying different reference frames. 
The procedure for searching the best possible values of a and /?, which would ensure convergence, is 
currently under investigation. The heuristics used for the tests (Frame change heuristics in section 
lOl) may certainly be improved taking into account the physical and geometrical features of the 
problem instance. 


Statistical results with change of frame 


Ariane launch 

Ariane flight 

Pegasus 

Number of cases 







- Total 

108 

234 

480 

- Easy (total execution time < 50 s) 

103 

120 

126 

- Diffi. (total execution time ^ 50 s) 


5 

114 

354 

Rate of success (%) 







- Total 

88.0 

54.7 

65.0 

- Before A 4 -continuation 

88.9 

88.0 

83.5 

- Average execution time (s) 

Easy 

Difh. 

Easy 

Diffi. 

Easy 

Diffi. 

- Total with change of frame 

2.43 

105.70 

11.38 

930.24 

32.6 

1122.2 

- Total with (a*,/3*) 

3.73 

13.33 

10.97 

61.45 

30.39 

50.01 

- In Ai-continuation 

0.22 

0.67 

0.62 

7.78 

0.52 

2.15 

- In A 2 -continuation 

2.63 

6.99 

4.79 

8.00 

6.56 

6.65 

- In Aa-continuation 

0 

0 

0 

0 

4.40 

12.23 

- In A 4 -continuation 

0.88 

5.67 

5.56 

45.67 

18.91 

28.98 

Average number of simulations 

Easy 

Diffi. 

Easy 

Diffi. 

Easy 

Diffi. 

- Total with (a*,/3*) 

123 

182 

122 

305 

246 

387 

- In the 1st continuation 

21 

25 

18 

54 

21 

32 

- In the 2nd continuation 

83 

125 

56 

55 

56 

52 

- In the 3rd continuation 

0 

0 

0 

0 

48 

150 

- In the 4th continuation 

20 

32 

47 

196 

122 

153 

Average times of change of frame 

1.4 

12.0 

0.0 

6.3 

0.2 

9.7 


Table 4: Statistic result with change of frame. 


Here we test the method with extreme conditions, whereas in real use, only small attitude 
maneuvers are used for Ariane-type launchers. Among the tested cases of Ariane launch (resp. 
Ariane flight), there are 86 cases with maneuver time t/ < 30 s and the average time to solve the 
problem is 6.85 s (resp. 3.76 s). For these small maneuvers, the continuation method provides 
rapidly an optimal control. 

Denote bill = |'0/ — V'ol the yaw change. We observe from the statistical results that, in 
case of Pegasus, the average execution time increases when Sijj increases. When (5'0 ^ 10°, the 
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average total execution time is 25.07 s, and when dij) ^ 20° the total average execution increases to 
110.30 s. Therefore, when applying to quasi-planar maneuvers {Sil) ^ 10°) for Pegasus like airborne 
launchers, the continuation procedure is also efficient. 

5.2 Comparison with the direct approach 

The proposed continuation method is much faster and more precise than the direct approaches. 
Furthermore, the trajectories obtained, even though they are sub-optimal, are physically feasible 
because the control remains smooth. On the other hand, the direct methods may produce oscil¬ 
latory control if chattering. In this section, we will use the direct method proposed in [15]. In 
the direct method, the midpoint rule (2nd order method) is used as discretization method and the 
number of time steps is set to be 100. Note that in the indirect method, we use Runge-Kutta 4th 
method which is more precise. 

Case without chattering We report on Figures andean optimal solution of (MTCP) ob¬ 
tained from the proposed continuation procedure. Here we consider the Ariane flight configuration 
in a planar case with 6*o = 38°, Of = 40° and Vg = 2000 m/s. Denote also the module of control as 
u, and thus ui = ucost/ and U 2 = usin^ with C defined as the angle of the control. 




Figure 2: State x{t) 


(continuation method). 


-R 0.5 


200 


10 20 
t: s 


30 



Figure 3: Control u{t) (continuation method). 


For this example, it takes about 8 s to solve the problem using the continuation procedure, 
while it takes about 40 s using the direct method initialized with an arbitrary constant commands. 
The optimal solution derived from the direct method is the same with the one reported on the 
Figures [D and E] except for some small oscillations around zero (namely for the variables ujx, 4’ 
and (/>). Moreover, the maneuver time tf = 30.03 s given from the direct method is about the same 
with the maneuver time tf = 30.07 s obtained from the continuation procedure. 
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From Figure 13] we see that the obtained optimal control has only two switches (at time 7.5 s 
and 22.7 s the control angle changes an angle of tt, and thus Ux changes from +1 to —1 at 7.5 s 
and changes from —1 to +1 at 22.7s). It is worth noting that for a real launcher, it is better to 
have a continuous control, since the acceleration of the control is generally limited. For this aim, it 
suffices to stop the A 4 -continuation before 1 to obtain a sub-optimal control. For example, we stop 
at 0.95 and get a control in Figure S) This control is feasible on a real launcher and the maneuver 
time is tf = 35.25 s, five seconds longer than the optimal one. 




Figure 4: Control u{t) (continuation method). 


Case with chattering In Figures |5| and | 6 | a sub-optimal solution obtained from the direct 
method is illustrated, while in Figures |7| and [Sj a sub-optimal solution is given by the continuation 
procedure (parameter A 4 stops at 0.994). Here we still consider the Ariane flight configuration in 
a planar case with 6 q = 30°, 9f = 40° and Uq = 3200 m/s. 



Figure 5: State x{t) and control u{t) (direct method). 



t: s 


200 



-200 -- 

O 50 

t: s 


Figure 6 : Control u{t) (direct method). 
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Figure 7: State x{t) and control u{t) (continuation method). 




Figure 8 : Control u{t) (continuation method). 


When using the continuation method (resp. direct method), the maneuver time obtained is 
tf = 61.04s (resp. tf = 58.64s) and the problem is solved within 41s (resp. within 182s). 
We observe that though the maneuver time derived from the direct method is less than with the 
continuation procedure, the sub-optimal control (see FigurelHl) given by the direct method oscillates 
much with a module less than 1 : this indicates that there is a singular arc in the “true” optimal 
trajectory, causing chattering at the junction with regular arcs. Such solution, with infinitely many 
control oscillations, cannot be used in the practice. 

From Figure [3 we observe that the sub-optimal control obtained from the continuation proce¬ 
dure is continuous and piecewise smooth. In fact, if we observe the time history of the switching 
function, we will see that the switching function tends to have more switches because of the chat¬ 
tering, and this cause the fail of the continuation method when A 4 > 0.995. 

In practice, the attitude maneuvers remain of limited magnitude for a launcher. The chattering 
phenomenon is thus unlikely to occur, excepted in very specific configurations, for example if the 
available nozzle deflection is very small. 


6 Conclusion 

In this paper, we have addressed the problem of the minimum time maneuver coupling attitude 
and trajectory dynamics for a launcher. The solution method is based on an indirect approach 
with a four-stage continuation procedure to solve the optimality conditions derived from the PMP. 
Starting from the explicit solution of a simplified problem of order zero, the successive continuations 
consists of retrieving the true attitude dynamics and terminal conditions. A regularization problem 
is used as an intermediate problem to overcome the problem caused by chattering. With this 
procedure, the solution to any problem instance can be obtained from scratch in a quite short time 
whatever the launcher data and the terminal conditions. The minimum time problem is fully solved 
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if possible. In case of a failure due to chattering, a very good sub-optimal solution is available 
using smooth control law suited to practical application. 

The method is exemplified on three representative launcher configurations, respectively for a 
vertical launch phase, an in-flight phase and an airborne launch phase. Statistical tests are run 
sweeping on the initial and final conditions. The statistical results show that the developed method 
is fast and robust. 

The current improvement axes deal with the choice of the coordinate frame that could be chosen 
from a physical and geometrical considerations. Moreover, state constraints on the angular velocity 
and the attack angle will be considered to provide safer optimal or sub-optimal trajectories for the 
launchers. 
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